#include "math.h"
#include <algorithm>
#include <string>
#include <sstream>
#include <fstream>
#include <iostream>
#include <vector>
#include <map>
#include <sys/time.h>
#include <unistd.h>
#include "Matrix.h"
#include <limits>

using namespace Numeric_lib;
using namespace std;

// возвращает текущее время в формате sec.msec
double cur_time() {
  struct timeval tv;
  gettimeofday(&tv, 0);
  return (double)tv.tv_sec + (double)tv.tv_usec*0.000001;
}

struct MICell {
    int m;
    int i;
};

int _nw_score(int reward, int penalty, int ogap, int egap,
              const string& seq1, const string& seq2,
              Matrix<MICell, 2>& matrix) {
    /*
    reward, penalty - веса за совпадение и несовпадение символов
    ogap, egap      - веса за открытие и продолжение разрыва
    seq1, seq2      - выравниваемые последовательности
    matrix          - выделенная структура в памяти для матрицы замен
    подробнее см. Biological Sequence analysis, Richard Durbin, p. 55
    */
    int i,j,s;
    int l1 = seq1.size();
    int l2 = seq2.size();
    MICell cell1, cell2;
    
    // initializing
    matrix[0][0].m = 0;
    matrix[0][0].i = 0;
    for (i=1; i<=l1; i++) {
        matrix[i][0].m = - ogap - (i - 1) * egap;
        matrix[i][0].i = - ogap - (i - 1) * egap;
    }
    for (i=1; i<=l2; i++) {
        matrix[0][i].m = - ogap - (i - 1) * egap;
        matrix[0][i].i = - ogap - (i - 1) * egap;
    }
    
    // calculate score
    for (i=1; i<=l1; i++) {
        for (j=1; j<=l2; j++) {
            if (seq1[i-1] == seq2[j-1]) {
                s = reward;
            }
            else {
                s = - penalty;
            }
            cell1 = matrix[i-1][j-1];
            matrix[i][j].m = s + max(cell1.m, cell1.i);
            cell1 = matrix[i-1][j];
            cell2 = matrix[i][j-1];
            matrix[i][j].i = max(max(cell1.m-ogap, cell1.i-egap), max(cell2.m-ogap, cell2.i-egap));
        }
    }
    return max(matrix[l1][l2].m, matrix[l1][l2].i);
}

int _nw_score_restricted(int reward, int penalty, int ogap, int egap,
              const string& seq1, const string& seq2,
              MICell* matrix) {
    /*
    reward, penalty - веса за совпадение и несовпадение символов
    ogap, egap      - веса за открытие и продолжение разрыва
    seq1, seq2      - выравниваемые последовательности
    matrix          - выделенная структура в памяти для матрицы замен
    подробнее см. Biological Sequence analysis, Richard Durbin, p. 55
    */
    int i,j,s,t1,t2;
    int l1 = seq1.size();
    int l2 = seq2.size();
    int k = 100;
    int MIN_INT = numeric_limits<int>::min();
    MICell cell;

    int max_i1, max_i2;
    
    // initializing
    matrix[0].m = 0;
    matrix[0].i = 0;
    for (i=1; i<=k; i++) {
        matrix[i*(l2+1)].m = - ogap - (i - 1) * egap;
        matrix[i*(l2+1)].i = - ogap - (i - 1) * egap;
        matrix[i].m = - ogap - (i - 1) * egap;
        matrix[i].i = - ogap - (i - 1) * egap;
    }

    // calculate score
    int tmp_index;
    for (i=1; i<=l1; i++) {
        t1 = max(1, i-k);
        t2 = min(l2, i+k);
        for (j=t1; j<=t2; j++) {
            if (seq1[i-1] == seq2[j-1]) {
                s = reward;
            }
            else {
                s = - penalty;
            }
            tmp_index = i*(l2+1)+j;
            cell = matrix[tmp_index-l2-2];
            matrix[tmp_index].m = s + max(cell.m, cell.i);
            
            if (j<i+k) {
                cell = matrix[tmp_index-l2];
                max_i1 = max(cell.m-ogap, cell.i-egap);
            } else {
                max_i1 = MIN_INT;
            }
            if (j>i-k) {
                cell = matrix[tmp_index-1];
                max_i2 = max(cell.m-ogap, cell.i-egap);
            } else {
                max_i2 = MIN_INT;
            }
            matrix[tmp_index].i = max(max_i1, max_i2);
        }
    }
    return max(matrix[l1*(l2+1)+l2].m, matrix[l1*(l2+1)+l2].i);
}

double nw_score(string seq1, string seq2, int reward, int penalty, int ogap, int egap) {
    int l1 = seq1.size();
    int l2 = seq2.size();
    Matrix<MICell, 2> matrix(l1+1, l2+1);
    
    return _nw_score(reward, penalty, ogap, egap, seq1, seq2, matrix);
}

double nw_score_mk(int reward, int penalty, int ogap, int egap,
                   string seq1, string seq2,
                   MICell* matrix,
                   int iter_num) {
    /*
    reward, penalty, ogap, egap - см. _nw_score()
    seq1, seq2 - выравниваемые последовательности
    matrix - выделенная структура в памяти для матрицы замен
    iter_num - число итераций в Монте-Карло
    */
    
    int i;
    int real_score;
    int rand_score;
    double scores_sum = 0;
    vector<int> scores(iter_num);

    real_score = _nw_score_restricted(reward, penalty, ogap, egap, seq1, seq2, matrix);
    // Monte-Carlo
    for (i=0; i<iter_num; i++) {
        random_shuffle(seq1.begin(), seq1.end());
        rand_score = _nw_score_restricted(reward, penalty, ogap, egap, seq1, seq2, matrix);
        scores_sum += rand_score;
        scores[i] = rand_score;
    }
    
    // calculating Z
    double ev = scores_sum / iter_num;
    double diff_sum = 0;
    for (i=0; i<iter_num; i++) { diff_sum += (scores[i] - ev)*(scores[i] - ev); }
    double sd = sqrt(diff_sum/iter_num);

    cout << "MK preview: " << real_score << "\t" << ev << "\t" << sd << "\n";
    
    if (sd > 0) {
        return (real_score - ev) / sd;
    } else {
        return 0.0;
    }
}
